library(haven)
library(tibble)

source("pests.r")
set.seed(9860)

mydata <- read_dta("replication data.dta")
mydata <- data.matrix(mydata)


# DIFFERENT LAG ORDERS
baseline.Par1 <- Parp(mydata[,"article_all"] ~ mydata[,"nat_cycle"] + mydata[,"newsp_fe1"]
 + mydata[,"newsp_fe2"] + mydata[,"newsp_fe3"] + mydata[,"newsp_fe4"] + mydata[,"newsp_fe5"]
 + mydata[,"month_fe1"] + mydata[,"month_fe2"] + mydata[,"month_fe3"] + mydata[,"month_fe4"]
 + mydata[,"month_fe5"] + mydata[,"month_fe6"] + mydata[,"month_fe7"]+ mydata[,"month_fe8"]
 + mydata[,"month_fe9"] + mydata[,"month_fe10"] + mydata[,"month_fe11"] + mydata[,"v15"]
 + mydata[,"lifts_restr_nat"] + mydata[,"sim_over"] + mydata[,"gov2"] + mydata[,"gov3"], p=1)
summary(baseline.Par1)

baseline.Par2 <- Parp(mydata[,"article_all"] ~ mydata[,"nat_cycle"] + mydata[,"newsp_fe1"]
 + mydata[,"newsp_fe2"] + mydata[,"newsp_fe3"] + mydata[,"newsp_fe4"] + mydata[,"newsp_fe5"]
 + mydata[,"month_fe1"] + mydata[,"month_fe2"] + mydata[,"month_fe3"] + mydata[,"month_fe4"]
 + mydata[,"month_fe5"] + mydata[,"month_fe6"] + mydata[,"month_fe7"]+ mydata[,"month_fe8"]
 + mydata[,"month_fe9"] + mydata[,"month_fe10"] + mydata[,"month_fe11"] + mydata[,"v15"]
 + mydata[,"lifts_restr_nat"] + mydata[,"sim_over"] + mydata[,"gov2"] + mydata[,"gov3"], p=2)
summary(baseline.Par2)

baseline.Par3 <- Parp(mydata[,"article_all"] ~ mydata[,"nat_cycle"] + mydata[,"newsp_fe1"]
 + mydata[,"newsp_fe2"] + mydata[,"newsp_fe3"] + mydata[,"newsp_fe4"] + mydata[,"newsp_fe5"]
 + mydata[,"month_fe1"] + mydata[,"month_fe2"] + mydata[,"month_fe3"] + mydata[,"month_fe4"]
 + mydata[,"month_fe5"] + mydata[,"month_fe6"] + mydata[,"month_fe7"]+ mydata[,"month_fe8"]
 + mydata[,"month_fe9"] + mydata[,"month_fe10"] + mydata[,"month_fe11"] + mydata[,"v15"]
 + mydata[,"lifts_restr_nat"] + mydata[,"sim_over"] + mydata[,"gov2"] + mydata[,"gov3"], p=3)
summary(baseline.Par3)
